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Abstract. We constrain the energy at which the transition from Galactic to extragalactic 
cosmic rays occurs by computing the anisotropy at Earth of cosmic rays emitted by Galactic 
sources. Since the diffusion approximation starts to loose its validity for E/Z > 10 16-17 eV, 
we propagate individual cosmic rays using Galactic magnetic field models and taking into 
account both their regular and turbulent components. The turbulent field is generated on 
a nested grid which allows spatial resolution down to fractions of a parsec. Assuming suffi- 
ciently frequent Galactic CR sources, the dipole amplitude computed for a mostly light or 
intermediate primary composition exceeds the dipole bounds measured by the Auger col- 
laboration around E ~ 10 18 eV. Therefore, a transition at the ankle or above would require 
a heavy composition or a rather extreme Galactic magnetic field with strength > 10 fiG. 
Moreover, the fast rising proton contribution suggested by KASCADE-Grande data between 
10 eV and 10 18 eV should be of extragalactic origin. In case heavy nuclei dominate the flux 
at E > 10 18 eV, the transition energy can be close to the ankle, if Galactic CRs are produced 
by sufficiently frequent transients as e.g. magnetars. 
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1 Introduction 

The question at which energy the transition from Galactic to extragalactic cosmic rays (CRs) 
takes place is one of the major unresolved issues of cosmic ray physics. Two promising 
possibilities are to associate the transition with one of the two evident features of the cosmic 
ray spectrum: The second knee around E ~ 5 x 10 17 eV or the ankle at E ~ 3 x 10 18 eV. 
Since the chemical composition of galactic and extragalactic CRs should differ in general, 
both because of propagation effects and of the different nature of their sources, the transition 
may be detected experimentally studying the chemical composition of CRs as function of 
energy. 

In the case of a transition around the second knee, Galactic CR sources such as e.g. 
supernova remnants would accelerate CRs up to the rigidity-dependent knee, which is close 
to 10 17 eV for iron. If the extragalactic CR flux dominating at higher energies would consist 
mainly of protons, the ankle could be explained as a dip in the extragalactic CR spectrum due 
to the pair-production losses of protons on cosmic microwave background (CMB) photons 
P+Tcmb — >• p+e + + e~ [1]. Below ~ 10 17 ~ 18 eV, the extragalactic CR flux may be suppressed 
because of CR propagation in extragalactic magnetic fields [2, 3]. On the other hand, the 
scenario of Ref. [4] would favour a transition at the ankle. The composition of the CR flux 
at high energies is the subject of current debate due to the facts that hadronic physics must 
be extrapolated from lower energies and that the complex experimental analyses for different 
experiments are not yet completely reconciled. The scenario of Ref. [1] is supported by the 
composition measurements of HiRes [5] and the first results of the Telescope Array [6], which 
are consistent with a light composition around the ankle and above. On the other hand, recent 
results from the Pierre Auger Observatory [7, 8] indicate a composition becoming heavier with 
increasing energy above the ankle, and the Yakutsk EAS array muon data suggests a non 
negligible fraction of heavy nuclei above ~ 10 19 eV [9]. Moreover, the measurements of the 
KASCADE-Grande [10] collaboration are consistent with a dominantly heavy composition up 
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to 10 eV. However, the KASCADE-Grande data indicate a fast rising proton contribution 
above 10 17 eV. 

Thus at present the experimental data on the CR composition do not allow us yet 
to determine the transition energy between Galactic and extragalactic CRs. In this paper 
we suggest to use instead experimental limits on the anisotropy of the arrival directions of 
UHECRs to constrain the maximal contribution of Galactic CRs at E > 10 18 eV. At energies 
below 10 17 eV, the diffusive propagation of Galactic cosmic rays and their resulting anisotropy 
at Earth was studied in details in Refs. [11, 12]. 

Since the propagation of CRs in the Galactic magnetic field (GMF) is not longer diffu- 
sive at E > 10 17 eV, we directly propagate UHECRs in the GMF using the numerical code 
developed in Refs. [13, 14]. We present also a way to generate the turbulent field on a nested 
grid without limitation on its spatial resolution. This method allows us to include mag- 
netic field fluctuations spanning the required large dynamical range of scales, from negligible 
compared to the CR Larmor radii up to 300 pc. As main result of this work we show that 
the existing limits on CR anisotropies strongly restrict the contribution of the CNO element 
group to the Galactic CR component above E > 1 EeV, while the contribution of iron is 
restricted above E > 3 EeV. 

Details of the method to generate turbulent magnetic fields are discussed in the Sec- 
tion 2. In Section 3, we review the GMF models used and discuss how the CR anisotropy 
is calculated. Results of numerical simulations are presented in the Sections 4 and 5 for 
anisotropies and the spectrum of UHECR. 



2 Modeling Turbulent Magnetic Fields 

We adopt in this section a convenient way to generate turbulent magnetic fields on nested 
grids which allows to include a large dynamic range of spatial scales contributing to the 
turbulence. 

A turbulent magnetic field B satisfies (B(r)} = and (B(r) 2 ^ = B 2 ms > 0. Let us 
denote k the modulus of wave vectors and a the spectral index of the field: a = 5/3, 3/2 
and 1 respectively for Kolmogorov, Kraichnan and Bohm spectra. The power spectrum of 
the field satisfies V(k) oc k~ a , and the amplitudes of its Fourier modes are \B(k)\ 2 oc k~ a ~ 2 . 
The spectral index of the turbulent Galactic magnetic field is poorly constrained. While 
a ~ 1 appears hardly plausible, both Kolmogorov and Kraichnan spectra could be allowed 
by the data. To study the dependence of our results on the spectral index, we present below 
computations for a = 5/3 and 3/2, as examples. Wave vector moduli satisfy 2-7r/L max < 
k = |k| < 2-7r/L m i n , where L m [ n and L max are respectively the minimal and the maximal 
variation scales in the turbulent field. In practice, L m i n corresponds to the damping scale of 
the field, which could be as low as an astronomical unit. We choose here L m i n = 1 AU. For 
a = 5/3 and 3/2, the value of L rrmi does not noticeably affect the results, because the larger 
a is, the more the energy is concentrated in the modes with large spatial variations. We take 
I max = 100 — 300 pc. The correlation length L c of the field, defined as in [15], is equal to 

j -^max Ot — — (-C'min/i-'max) (o -\\ 

c ~~ ~~ 9 1 _ (T ■ IT V* 1 ! ' ^ ' 

As discussed in Refs. [14, 16], there are two main numerical methods to generate tur- 
bulent magnetic fields. First, they can be generated as a superposition of plane waves as 
in Ref. [17] and computed in any point of the space. Second, values of the field can be 
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pre-computed with the Fast Fourier Transform on a three dimensional cubic grid, which is 
periodically repeated in space. The value of the field can be extrapolated to any position 
from these values. Computing the individual trajectories of millions of cosmic rays with 
rigidities as low as E/Z ~ 3 x 10 16 eV is achievable within reasonable computing times only 
with the second method. The number of vertices on such cubic grids is AA 3 . J\f ~ 256 — 512 is 
typically the limit above which the grid cannot be loaded in a 2 gigabyte RAM memory. The 
ratio L max /L m ; n is limited by N/2, when L max equals the size of the cubic box. Moreover, 
we take L mSuX / X m ; n to be smaller than N/2, by at least a factor of a few. This ensures that 
the modes with the largest spatial variations ~ L max have a few oscillations within the box 
size. Otherwise, the generated turbulent field can be highly anisotropic. Cosmic rays which 
diffuse in turbulent magnetic fields are mostly sensitive to modes with wave numbers k close 
to ~ 27t/Vl, where tl is their Larmor radius. For E/Z = 10 18 /26eV and a field of strength 
6/iG, it is equal to tl — 7pc. In the numerical simulations, one can disregard modes with 
2ir/k <C ri, because they have a negligible influence on the particle trajectories. On the con- 
trary, modes with 2ir/k € [r^, L max ] which isotropize cosmic rays in a non trivial way have to 
be taken into account. Therefore, instead of using 2ir/k £ [L m i n , L max ] = [1 AU, 100 — 300 pc] , 
we truncate the minimal scale of spatial variations for the generated field and restrict our- 
selves to [L' min ,L max ] with L' mill sufficiently small compared to However, L ma _ x /L' min is 
still too large to fit in one magnetic field grid of reasonable size. To solve this issue, we use 
the method of nested grids explained in the following. 

Let us assume that B(r) is the sum of N + 1 components: B(r) = X^o-^*( r ) (^ or 
j ^ i, (Bi(r) ■ Bj(r)) = 0). In practice, TV = 2 is sufficient for this work. Bq{t) contains 
all Fourier modes with 2n/k € [L mm ,L min ], and the fields Bi(r) (1 < i < N) respectively 
contain the modes with 2-n/k € [Lj,Lj + i], where L\ = L' min and Ljy + i = L max . The ratios 
Lj/Lj+i are all chosen to be smaller than N/2 by a factor of a few. 

The root mean square (rms) strength of the total field, -B rms , satisfies [18] 

p2ir/L min 

B? ms oc / dkV(k) . (2.2) 

max 



Therefore, for a ^ 1, £ r 2 ms cx (L m ~ x - L^ 1 ). For < i < N, the energy density present in 
Bi(r) is proportional to 



max mm 

Lmax ^min 



which yields the rms amplitude of Bi(r), B TmSt i. The turbulent field -B tur b generated for the 
computations is equal to the sum of the N components Bi(r) with i = 1,...,N, -B tur b = 
^2^=1 Bi(r). -B tur b is equal to the total turbulent field B(r) after subtracting the modes 
with spatial variation scales smaller than i min <C r^. Each Bi(r) is generated on a cubic 
grid of lateral size NLi/2. Each grid is periodically repeated in physical space. The Bi(r) 
with large i contain the modes with large spatial variation scales and the Bi(r) with small 
i, the modes with small variation scales. In any space point, the magnetic field from the 
large (respectively small) resolution grid is evaluated as the 8-point linear interpolation of 
the values on vertices of the large (respectively small) scale resolution grid. 

We have verified that we recover with this code the results found by the earlier studies 
of Refs. [16, 20]. As an example, we present in the appendix our computations of the CR 
diffusion coefficient for pure magnetic turbulence, as well as the parallel and perpendicular 
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diffusion coefficients for turbulence superimposed to a regular field. Our results are found to 
be in very good agreement with those of these previous studies. 

To summarize, we take in the following L m in = 1 AU for the normalisation of the 
turbulent field strength, so that the rms for the total field with modes satisfying 2n/k £ 
[Lmin = 1 AU, L max ] would be B rms . In all computations for CRs with rigidities E/Z < 1 EeV, 
we take N = 2 and set L' min , the actual minimal scale of fluctuations in the generated field, 
to 1 pc. The intermediate scale L2 between the N = 2 grids is 20 pc. In practice, for large 
rigidities E/Z > lEeV, tl > 180 pc and we can drop the smaller scale grid and only use 
the larger one : N = 1 and L' min = 20 pc. We use the standard Runge Kutta method with 
adaptative step size of Ref. [19]. Removing the smaller scale grid increases the step size of 
the integrator and allows us to reduce computing time for particles with E/Z > lEeV. 



3 Galactic Magnetic Field models and Method to compute the Anisotropy 

The Galactic magnetic field (GMF) can be regarded as the sum of a regular component (large 
scale variations) and a turbulent component (small scale variations). 

We described in the previous Section a method to generate numerically the turbulent 
component. The spatial profile of the rms strength of the turbulent field, B rras (r, z), is poorly 
constrained. Therefore, we use two different types of profiles as examples. First, we take a 
model with an exponentially decaying field strength in the Galactic halo [21]. We will refer 
to it as the "Profile 1" : 

B rms (r,*)=B(r)expf-Mj , (3.1) 

where r is the Galactocentric radius and z the distance to the Galactic plane. The parameter 
zq denotes the scale height of the random field into the z— direction. We will take zq = 
(2 — 8) kpc in this work. The radial profile B(r) is equal to 

/ B o«P(M) , if r<3kpc (bulge) 

B(r) = l^ex P (^-te.),if r >3kpc < 3 ' 2 > 

where Bq is defined as the value of B rms close to the Sun. 

Second, we also consider a constant rms strength within a box of size r < 20 kpc and 
\z\<z ("Profile 2"): 

R (V ~\ - / B > if r ^ 20k P C and M < Z (O o\ 

[0 , if r > 20 kpc or \z\ > zq 

Although this profile is very likely less realistic than the previous one, we test it because it 
corresponds to the profile used in the usual "leaky-box approximation" . 

The global geometry of the regular GMF is still poorly known. The Faraday rotation 
measures (RM) for extragalactic sources suggest that it is made of at least two different 
components, in the disk and in the halo, with different geometries [22, 23]. The field in the 
disk is believed to be symmetric with respect to the Galactic plane, while the field in the 
halo is believed to be antisymmetric [22, 23]. The RM at high latitudes show the existence of 
a toroidal field in the halo, on each side of the Galactic plane. This field is counter clockwise 
in the Northern halo and clockwise in the Southern halo, as seen from the Galactic North 
pole. Several analytical models have been proposed to describe the regular GMF. As shown 
in Refs. [22, 24], presently no theoretical GMF model can fit all experimental data. However, 
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Ref. [23] presents the two first models that fit reasonably well all extragalactic RM data in 
most regions of the sky, which represents a significant improvement of our knowledge of the 
GMF. Since it is impossible to distinguish between bisymmetric (BSS) and antisymmetric 
(ASS) geometries for the disk field, the authors of Ref. [23] propose two different benchmark 
models for the regular GMF. Below, we will refer to them as the "PTKN-BSS" and "PTKN- 
ASS" models. They contain disk and toroidal contributions. Let us use Galactocentric 
cylindrical coordinates (r,9,z), where r = (rc 2 + y 2 ) 1 ^ 2 , and Cartesian coordinates x, y and 
z. The Earth is assumed to be at (x = 0, y = r & = 8.5 kpc, z = 0), where 9 is set to zero 
at the position of the Earth, and increases clockwise, as seen from the Galactic North pole. 
The components in cylindrical coordinates of the disk field strength, B r and Bq, are defined 



as 



B T = B (r, 6, z) ship, 
Bg = B (r, 8, z) cosp , 



(3.4) 



with p = —5° and p = —6° respectively for the ASS and BSS models. For the ASS disk field, 



cos 



B(r,9,z) = b(r 
while, for the BSS disk field, 

B (r, 9, z) = b (r) cos 



1 



hi 



+ 



cxp 



tanp 



In ( — ] + 

re 



\z\ 



•exp 



(3.5) 



(3.6) 



where (/>=!/ tanp ■ ln(l + d/r & ) — ir/2 with zq = 1.0 kpc, d = —0.6 kpc and 



b(r) 



5 . ok pc Q co S0 forr<5.0kpc 
2 -°^ G 7^b; for r> 5.0 kpc 



The halo field components -Btx and -&ry are defined as 

i?Tx = — -St sgn (z) cos 6, 
Btj = -Bt sgn (z) sin 6 , 



where 



exp ' 



TO 



1 + 



\z\-h T 

Wt 



(3.7) 
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with -Btoi r T0 ; h'Y chosen as in Table 1. For \z\ < h^, wt = 0.25 kpc, and for \z\ > Jit, 
w<t = 0.4 kpc. The strength of the halo field decays towards the Galactic center, for r < rxo- 

For most of the following computations, we use the PTKN-BSS model as an example. 
We test the dependence of our results on the regular GMF by also using the PTKN-ASS 
model, the "ASS+RING" model of Ref. [25] (which we will refer to as "Sun08" in the follow- 
ing), and the Prouza and Smida (PS) model [26, 27] with the parameters given in Ref. [13]. 

The Galactic center [28-30], some types of supernovae [31], magnetars [32-35] or GRBs [36- 
42] have been discussed as potential Galactic sources able to accelerate CRs to ultra-high 
energies. The spatial extension of the region containing Galactic CR sources is better con- 
strained than the GMF parameters. Sources are expected to be distributed in the Galactic 
disk, within ~ ±(200-500) pc from the Galactic plane z = [43, 44]. The Galactocentric 
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ASS, z > 


ASS, z < 


BSS, z > 


BSS, z < 






2fiG 


4/xG 


4^G 




6 kpc 


6 kpc 


6 kpc 


5 kpc 




1.3 kpc 


1.3 kpc 


1.5 kpc 


1.5 kpc 



Table 1. Values for the Northern [z > 0) and Southern (z < 0) halo parameters Bto, ^to and /it, 
in the ASS and BSS versions of the PTKN model. 




Figure 1. Trajectories of iron anti- nuclei backtraccd from the Earth. Left panel: Energies equal 
to (1, 4, 8) x 10 18 eV; Right panel: Energies equal to (1, 4, 10) x 10 19 eV. For details on the Galactic 
magnetic field model, see text. 

radius r up to which the source region extends is less constrained. Since there should not be 
a significant number of sources with r > 20 kpc [43], we will take in the following, for most 
cases, r = 20 kpc as the limit of the source region. 

At sufficiently low rigidities, the Larmor radius of cosmic rays is smaller than the co- 
herence length of the turbulent GMF. Previous studies that predicted the amplitude of the 
cosmic ray anisotropy at Earth assumed CRs are diffusive. While the diffusion approxima- 
tion is justified for rigidities smaller than E/Z ~ 10 17 eV, it starts to fail in the rigidity 
range investigated in this work: E/Z > (10 18 /26) eV. Between these rigidities, one typically 
expects a transition from the diffusive regime to the ballistic regime for CR propagation. The 
transition does not happen abruptly at a given rigidity, which leads to non-trivial modes of 
CR propagation. This can have a non-trivial impact on the anisotropy of Galactic CRs at 
Earth. Therefore, we propagate in this work individual cosmic rays in models of the GMF. 

Figure 1 shows trajectories of anti-iron nuclei with energies 10 18 eV < E < 10 20 eV 
backtraced in one GMF model from the Earth, located at {x = 0, y = 8.5 kpc, z = 0). It 
shows the variety of CR propagation types in the transition from "purely" diffusive (here at 
10 18 eV) to "purely" ballistic (here above > (2 - 4) x 10 19 eV). This energy range is shifted 
when the magnetic field parameters are changed. The regular GMF used for these plots is the 
PTKN-BSS model. The turbulent component has a Kolmogorov spectrum with L m [ n = 1 AU, 
L max = 200 pc, the profile 1 with zq = 2 kpc, and a strength set to B IU1S = 4//G. 

The left panel of Figure 1 display the trajectories of 1, 4 and 8 xl0 18 eV iron anti- 
nuclei. Values of spatial coordinates on the axes are given in kilo-parsecs. For these GMF 
parameters, the Larmor radius of the 10 18 eV nuclei is smaller than the correlation length 
L c ~ 40 pc of the turbulent component. The trajectory of this cosmic ray resembles a random 
walk, see the red line. The green (4 x 10 18 eV) and blue (8 x 10 18 eV) trajectories, respectively, 
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correspond to diffusion in the regimes when tl — L c and tl > L c . The trajectory of the 
8 x 10 18 eV iron anti-nucleus is still confined in the Galactic plane for an extended time. On 
the right panel of Figure 1, one can see that this anti- nucleus goes back and forth in the disk. 
It propagates especially along the regular field lines which are locally approximately oriented 
along the x axis. 

The right panel of Figure 1 shows the trajectories of 1, 4 and 10 x 10 19 eV iron anti-nuclei. 
At 10 19 eV (red line), the CR is still strongly deflected before escaping the Galaxy. If one 
sums up all deflections along its trajectory, it exceeds 360°. This iron anti- nucleus is weakly 
deflected over distances up to ~ 1 kpc. It is strongly deflected only locally, when it reaches 
regions with stronger turbulent magnetic field fluctuations. At 10 20 eV, the trajectories are 
fully ballistic, see the blue line. At such energies, one expects that iron nuclei suffer deflections 
of the order of ~ 20° — 40° before escaping the Galaxy. The 4 x 10 19 eV anti-nuclei are not 
diffusive any more. However, they still experience large deflections. For instance, there is a 
big wiggle on the 40EeV particle trajectory, see the green line in the right panel. 

In principle, the best way to compute the anisotropy of Galactic CR at Earth is to 
use forward tracking. One should inject cosmic rays in the Galaxy at the source locations 
and only record the momenta of CR which cross a sphere around the Earth. The radius 
of this sphere should be small compared to the CR Larmor radius. This is, however, not 
feasible within reasonable computing times for the lowest rigidities we study. Therefore, we 
use a method first proposed in Ref. [45], and reused in more recent works such as Ref. [46]. 
It consists in backtracking anti-particles with random initial momenta from the Earth to 
outside the Galaxy, and to record for each one the total path length in the source region. This 
corresponds to assuming a continuous and homogeneous source distribution inside the source 
region. Except for the Galactic center, potential Galactic CR sources should be transient. 
The method we use here assumes the existence of sufficiently frequent transient sources in 
the Galactic plane, so that the continuous source distribution hypothesis is fulfilled. For rare 
transient sources whose periodicity start to be comparable to the CR confinement time in the 
Galaxy, the current anisotropy may strongly differ from the average anisotropy [47] . Ref. [47] 
shows that if GRBs were to be the sources of Galactic CRs up to the ankle, strong variations 
of the Galactic CR flux and anisotropy should be expected on time scales of a hundred 
Myr. This would mean that CR anisotropy limits may be compatible with Pierre Auger 
measurements if we live in atypical times, when it drops below ~ 2% at EeV energies [48]. 
Since one can hardly go beyond this statement for rare transients such as GRBs, we restrict 
our work to the case of sufficiently frequent transients. We will discuss below in more detail 
the domain of validity of this approximation, see Section 5. 

We count the length of particle trajectories contained in the source region without any 
weighting depending on the position inside the region. This corresponds to assuming that 
the sources are homogeneously distributed inside — 200 pc < z < 200 pc and r < 20 kpc. 
In practice, taking a more realistic source distribution with, for example, a source density 
decreasing with \z\, would only increase the Galactic CR anisotropy at Earth. The values 
presented below can be regarded as lower limits. 

Finally, we compute the amplitude of the dipole of CR anisotropies and compare it to 
the upper limits presented by the Pierre Auger Collaboration, for energies E > 1 EeV [48] . 
To do so, we associate to each of the N backtracked cosmic rays a vector v(9,<f>) whose 
direction on the sky (6,(j>) corresponds to the initial CR direction at Earth. Its length \v\ 
corresponds to the trajectory path length in the source region. The dipole direction is given 
by the sum of all vectors J2i v i- I n spherical coordinates (r, 6, (j>), v = L (1 + V cos 6) u r for 



-7- 



a dipole of amplitude T> and direction u z , with L = jj \vi\ being the average length of 
vectors and 9 the angle between u r and u z . Since the vectors v are isotropically distributed, 



i=l 

If higher order multipoles are present, Eq. (3.9) is unchanged because only the dipole has a 
non-zero contribution to the sum ^ t>j. 

4 Anisotropy of Galactic Cosmic Rays predicted at Earth 

In this Section, we compute the anisotropy at Earth of cosmic rays emitted by sources 
distributed in the Galactic plane, and compare it to the upper limits on the anisotropy as 
measured by the Pierre Auger Collaboration [48]. Within a given Galactic magnetic field 
model, if the predicted anisotropy of Galactic cosmic rays exceeds these limits at a given 
energy, a sufficiently large contribution of extragalactic CRs is required above that energy. 
Extragalactic CRs of energies below the GZK threshold around ~ 4 x 10 19 eV [49, 50] tend 
to have anisotropics at the percent level or below because either a large number of sources at 
cosmological distances can contribute to the flux or, in case of relatively strong deflections 
in extragalactic magnetic fields, such deflections tend to wash out anisotropics over the long 
path lengths propagated over the typical energy loss distance. Resulting anisotropies of 
extragalactic cosmic rays can, therefore, easily be below the current upper limits. If sources 
at cosmological distances dominate the extragalactic CR flux, the motion of the Sun with 
respect to the CMB frame would induce an anisotropy of the extragalactic flux of ~ 0.6% [51], 
due to the Compton-Getting effect [52]. 

In subsection 4.1, we present our predictions for the Galactic CR dipolar anisotropy and 
discuss its implications on the energy at which the transition from Galactic to extragalactic 
CR should occur. Figs. 2-5 show computations for this anisotropy. Assuming in a first 
approximation that the extragalactic flux is isotropic, one can deduce from these figures 
the maximum contribution of Galactic CR sources to the total flux. For instance, at the 
transition from Galactic to extragalatic CRs, half of the CR flux is of Galactic origin and 
half of extragalactic origin. Therefore, the exact transition energy must be in the energy 
range where half of the dipole amplitude of Figs. 2-5 does not overshoot the experimental 
upper limits on it. In 4.2, we study in more detail the anisotropy predicted by Galactic 
protons below the ankle and its dependence on the GMF parameters. 

4.1 Dipole Amplitude and predicted Energy of the Transition 

We assume that the sources of Galactic CRs are located in a cylinder along the Galactic 
disc with height z mSjX and radius r max . Varying the extension of the source region from 
r max = 15kpc to r max = 20kpc, we verified that for rigidities E/Z < 3 x 10 18 eV, the dipole 
amplitude does not change by more than 3%. For higher rigidities, the difference rarely 
exceed (5 — 10)%. In practice, the conclusions below will not change for r < 15kpc or 
r < 20 kpc. In the following we assume r < 20 kpc. 




and the dipole strength equals 




(3.9) 
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Figure 2. Left panel: Predicted amplitude of the dipole as measured at Earth versus energy, for 
different primaries (p, He, C, Si, Fe) emitted by Galactic sources distributed in the region — 200 pc < 
z < 200 pc and r < 20kpc. The dashed blue line represents the 99% C.L. upper limit on the dipole 
amplitude in right ascension as measured by the Pierre Auger Observatory [48] . Blue points represents 
the Pierre Auger measurements of the dipole in right ascension with the "East- West method" for the 
1 — 2 EeV bin and with the "Rayleigh analysis" for the three other energy bins, according to Fig. 5 
of Ref. [48]. The PTKN-BSS model was assumed for the regular GMF. The turbulent component is 
assumed to have a strength Bq = 4/xG, profile 1, and zq = 2kpc for its extension into the halo, with 
limiting length scales L min = 1 AU and i max = 200 pc; Right panel: Same as for the left panel, but 
for profile 2. 

Larger extensions of the source region along the z— direction would reduce the predicted 
anisotropy. We verified that, as long as the source region is less extended than \z\ < 500 pc, 
no strong modification of the results and conclusions would arise. We assume \z\ = 200 pc 
for most of the figures below. 

For each of the following plots we backtrack 10 4 particles. We find that this induces 
an error < ±3% on the predictions of the dipole amplitude. Reducing further this error to 
< ±1% could in principle be achieved by backtracking ~ 10 times more CRs, but it is in 
practice impossible due to computing time reasons. 

Figure 2 presents the simulated predictions for the dipolar amplitude at Earth for several 
different Galactic CR primaries: protons, helium, carbon (representative for the CNO group), 
silicon and iron. CR sources are assumed to be located in the thin disk with — 200 pc < 
z < 200 pc and r < 20kpc. The PTKN-BSS model is assumed for the regular GMF, and 
the turbulent field parameters are taken to be Bq = 4/zG, zq = 2kpc, L m i n = 1 AU and 
L max = 200 pc. In the left panel, we use profile 1 given in Eq. (3.1) for the turbulent field, 
whereas in the right panel we take profile 2 from Eq. (3.3). We include in these figures 
measurements of, and 99% C.L. upper limits on, the dipole amplitude in right ascension 
from the Pierre Auger Collaboration, as indicated in the captions. We assume that the 
true dipole vector does not lie in the equatorial plane, so that these upper limits in right 
ascension do not overconstrain significantly the total amplitude. In Fig. 2, the predictions 
for the dipole amplitude for silicon and iron primaries may appear to exceed the Auger upper 
limits in the whole energy range we consider. However, since error bars on our computations 
are ~ ±3%, some parts of these silicon and iron lines are compatible with the Auger limits at 
low energies. For instance, one cannot exclude the dipole amplitude for 1-fewEeV Galactic 
iron to be below the Auger limits. In the next two figures, 3 and 4, the red dashed lines 
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represent the lower error bars on the red shaded bands. These bands represent the allowed 
ranges for the dipole amplitude with Galactic iron primaries, for different turbulent GMF 
parameters. In the energy ranges where this dashed line falls below Auger upper limits, 
a pure flux of Galactic iron CRs cannot be excluded on anisotropy grounds. In practice, 
we expect that with ten times larger statistics, predicted dipole amplitudes at the lowest 
rigidities (E/Z < fewEeV/26), which are already smaller than a few percent, would fall 
further below the Auger limits. We have 10 sets of 10 4 CRs for different turbulent GMF 
parameters. Their individual dipole directions look random at the lowest rigidities and if one 
adds up these ten sets, the resulting anisotropy for iron at E = 10 18 eV falls to ~ 0.6%, which 
is below the Auger limit. This gives an idea of what should be expected for ten times larger 
statistics. Therefore, below the ankle, both silicon and iron of Galactic origin are compatible 
with Pierre Auger Observatory limits. For the set of GMF parameters assumed here, in case 
of a predominantly heavy composition below the ankle and sufficiently frequent transient 
sources, CRs may still be of Galactic origin up to the ankle. Iron nuclei of Galactic origin 
up to ~ 10 EeV cannot currently be ruled out from the point of view of the CR anisotropy. 

Depending on the composition at E ~ 10 18 eV, this has an important implication for the 
transition energy between Galactic and extragalactic CRs: For the set of GMF parameters 
assumed here, if the CR primary composition is predominantly light (p, He) or intermediate 
(C, N, O) at these energies, the predicted anisotropy at Earth would be larger than the 99% 
C.L. upper limits from the Pierre Auger experiment if these nuclei were of Galactic origin, 
as seen in Fig. 2. This implies that if the composition at E ~ 10 18 eV is measured to be 
light or intermediate, scenarios in which the transition from Galactic to extragalactic CRs 
occurs at the ankle are strongly disfavoured, at least for a wide range of GMF parameters. 
We investigate below the ranges of parameters for which this conclusion would be valid. 

Figure 2 also shows that the conclusions above do not strongly depend on the turbulent 
field profile. For profile 2 (right panel), the predicted dipolar anisotropy grows slightly more 
slowly with energy than for profile 1 (left panel). This is expected because for profile 1 
the gradient of the turbulent field tends to drive CRs towards larger z in the Galactic halo 
slightly faster than for the constant field of profile 2. Since predicted anisotropies are not 
very different for the two profiles, we will mostly focus on profile 1 in the following. 

We have also tested the dependence of these results on the regular GMF model. For 
rigidities E/Z > 3 EeV, the dipole amplitude and direction depends on the regular GMF 
model. However, the change of the dipole amplitude is too small to affect significantly our 
findings. Moreover, the PTKN-BSS model which we use in all the following figures is one of 
the models with the lowest dipole amplitudes among those tested. 

Figures 3 and 4 present the dependence of the previous results on the other turbulent 
GMF parameters, which are poorly constrained: the index a of the fluctuation spectrum 
(Fig. 3 - left panel), the maximal length of field fluctuations L max (Fig. 3, right panel), the 
field strength normalization Bq (Fig. 4, left panel) and the scale height zq (Fig. 4 - right 
panel). The shaded areas of the filled cures represent the relative change of the results for p, 
C and Fe primaries when varying separately the four above parameters and keeping all other 
parameters at the values in Fig. 2. The red dashed line is computed as the lower boundary 
of shaded areas for iron primaries minus 3% from our statistical uncertainties. 

At low rigidities E/Z < 4 EeV, the dipole amplitude grows with E/Z as expected. At 
larger rigidities, CRs start to enter the ballistic regime and higher order multipoles start 
to make a significant contribution to the total Galactic CR anisotropy at Earth. As seen 
in Figs. 3 and 4, the dipole amplitude may then become smaller and/or vary with E/Z. 
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Figure 3. Dependence of the predicted dipole amplitude on the turbulent field spectral index and on 
its maximum spatial variation scale. For comparison, Pierre Auger data [48] are shown in blue as in 
Fig. 2. Left panel: Shaded area for a G [3/2,5/3] (from Kraichnan to Kolmogorov); Right panel: 
Shaded area for L max G [100 pc, 300 pc]. Red dashed lines for the lower error bars on the iron filled 
curve. For each plot, the values for all other parameters are set to those used in Fig. 2. 




Figure 4. Dependence of the predicted dipole amplitude on the turbulent field strength at Earth 
and on its extension in the Galactic halo. For comparison, Pierre Auger data [48] are shown in 
blue as in Fig. 2. Left panel: Shaded area for Bq G [2/iG,8^G]; Right panel: Shaded area for 
zq G [2kpc, 8kpc]. Red dashed lines for the lower error bars on the iron filled curve. For each plot, 
the values for all other parameters are set to those used in Fig. 2. 

Thus a decrease of the dipole amplitude at high energies does not necessarily imply that the 
distribution of CR arrival directions at Earth becomes more isotropic. 

The widths of the filled curves indicate that results mostly vary with the turbulent GMF 
strength and its maximum spatial variation scales. Results are less sensitive to the spectral 
index of the field and no strong difference in the dipole amplitude is found between the fields 
with Kolmogorov and Kraichnan spectra. The amplitudes at E/Z > 3EeV/26 are slightly 
larger for the Kraichnan spectrum because less power is present in the large length scale 
modes, which are relevant at such rigidities, than for the Kolmogorov spectrum. Results are 
only marginally sensitive to the extension zq of the turbulent field into the halo, see right 
panel of Fig. 4. This is due to the fact that CRs which escape the source region and propagate 
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Figure 5. Predicted amplitude of the dipole versus the turbulent Galactic magnetic held strength, 
for E/Z = 10 18 eV cosmic rays emitted by Galactic sources distributed in the region with r < 20kpc, 
and — 200 pc < z < 200 pc or — 500 pc < z < 500 pc, respectively, as indicated. Prohles 1 and 2 (see 
subsection 3) for the turbulent Galactic magnetic held prohle along z, as indicated. For the regular 
GMF the PTKN-BSS model is assumed. Shaded or delimited areas correspond to zq varying in the 
range 1 — 8kpc. For the turbulent component a Kolmogorov (left panel) or Kraichnan (right panel) 
spectrum with L m - m = 1 AU and L max = 200 pc is assumed. 



to large z in the halo rarely come back to the source region. 

Results are mostly sensitive to L max and Bq. The maximal length scale of the turbulence 
L max determines up to which energy CRs still scatter on the turbulent field inhomogeneities. 
For larger L max , CRs can be diffusive up to larger energies, which therefore reduces their 
expected anisotropy at Earth. One can see in the right panel of Figure 3 that for L max = 
300 pc and Bq = 4/iG, the dipole amplitude below E ~ 15EeV for a pure iron composition 
may be compatible with the current 99% C.L. upper limits from the Pierre Auger Observatory. 

The dependence of our results on the turbulent field strength Bq is very strong, see the 
left panel of Fig. 4. The upper parts of the shaded areas correspond to Bq = 2 fiG and the 
lower to So = 8/iG. For Bq = 8/iG, the anisotropy at Earth of iron primaries is a priori 
compatible with the Pierre Auger upper limits up to E ~ 20 EeV. For Bq = 2 /iG, the dipole 
amplitude starts to overshoot the Pierre Auger limits around E ~ 3 x 10 18 eV, while for 
So = 4/iG, the amplitude starts to exceed the Pierre Auger limits around 10 EeV. For all 
cases, a light or intermediate composition at E ~ 10 18 eV would exceed the Pierre Auger 
upper limits and rule out the ankle as the transition from Galactic to extragalactic cosmic 
rays. 

4.2 Dipole Amplitude at E/Z = 10 18 eV 

In this section we demonstrate that for any reasonable combination of turbulent GMF param- 
eters, at E ~ 10 18 eV the dipolar anisotropy predicted by light primaries of Galactic origin 
is always larger than the observational limit from the Pierre Auger experiment. Therefore, 
having a reliable composition measurements at such energies is crucial for knowing if the an- 
kle can or cannot be the signature of the transition from Galactic to extragalactic CRs. Low 
energy extensions such as HEAT [53] and AMIGA [54] can solve this important question. 

Figure 5 shows how the strength predicted at Earth of the dipole amplitude of 1 EeV 
protons from Galactic sources depends on the turbulent field rms strength Bq = 2 — 8 /iG 
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Figure 6. Histograms of fractions of cosmic rays spending a given time in the source region (—200 pc < 
z < 200 pc and r < 20kpc), for rigidities E/Z = (1, 1.78, 3.16) x 10 18 eV/26 (left panel), and 
E/Z = (5.62, 10, 17.8) x 10 18 eV/26 (right panel). For the regular GMF the PTKN-BSS model is 
assumed. For the turbulent component a Kolmogorov spectrum with L m in = 1 AU and L max = 200 pc, 
strength Bq = 4 /iG and scale height zn = 2 kpc is assumed. 

when the scale height zq is allowed to vary in the range 1 to 8 kpc (shaded areas). Both 
turbulent field profiles are tested, and a ±500 pc width source region is also tested. Both for 
Kolmogorov (left panel) and Kraichnan spectra (right panel) , predicted dipole amplitudes are 
above 10%, considerably higher than the ~ 2% upper limit from the Pierre Auger experiment 
at such energies. Therefore, the ankle cannot be the signature of the transition from Galactic 
to extragalactic CRs, if the contribution of Galactic protons to the CR flux at 1 EeV is larger 
than ~ 20%. Such a scenario would be consistent with the dip model [1]. A CR flux at 1 EeV 
consisting mostly of light nuclei cannot be of Galactic origin, except in the very unlikely case 
of B » 10 fiG. 

5 Energy Spectrum of Galactic Cosmic Rays and Sources contributing at 
Earth 

Figure 6 presents histograms of the relative fraction of CRs spending a certain time in the 
source region, for rigidities ranging from lEeV/26 to 17.8EeV/26. We use the PTKN-BSS 
model for the regular GMF component. The turbulent GMF strength is set to Bq = 4 fiG, its 
extension in the halo to zq = 2 kpc. We take a Kolmogorov spectrum with maximum spatial 
variation scale L max = 200 pc. 

With such parameters, the average time spent in the source region for cosmic rays with 
E/Z = lEeV/26 is ~ 0.5 Myr. The CR escape times from the magnetized region of the 
Galaxy defined as —10 kpc < z < 10 kpc and r < 20 kpc, are found to be ~ 5 times larger 
than the times spent in the source region for the turbulent GMF profile 1. CRs which escape 
the Galactic thin disk containing CR sources still stay a non-negligible amount of time in the 
halo compared to the time spent in the source region. 

For 1 EeV iron nuclei (red curve in Fig. 6 - left panel) , only a few percent of CRs stay 
more than 1 Myr in the source region. This implies that rare transient sources such as gamma 
ray bursts (GRBs) are very unlikely to be sources of Galactic CRs in the sub-ankle region, 
even if heavy nuclei were able to escape such sources. 
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Figure 7. Relative contributions per volume to the Galactic CR flux observed at Earth: Fraction 
of all particles backtraced from the Earth which cross (200 pc) 3 cubes located in the source region. 
Earth is located at (x,y) = (0,8.5 kpc). Galactic plane in the plane of the panels. Upper left 
panel: For rigidity E/Z = 10 18 cV/26; Upper right panel: E/Z = 3 x 10 18 eV/26; Lower left 
panel: E/Z = 10 19 eV/26; Lower right panel: E/Z = 3 x 10 19 eV/26. Same Galactic magnetic 
field parameters as in Fig. 6. 

For lOEeV iron nuclei, the average time spent in the source region is ~ 0.06 Myr, 
which is nearly ten times smaller than at 1 EeV. We found that the average time spent in 
the source region is approximately proportional to 1/E -or slightly softer- in the rigidity 
range E/Z G [lEeV/26, 20EeV/26]. This implies that for a source injection spectrum 
proportional to E a , the spectrum at Earth approximately goes as E 01-1 . Interestingly, this 
result is compatible with predictions from diffusion in the regime where the Larmor radius is 
larger than the coherence length of the turbulent field: The distance traveled by CRs within 
a given amount of time is proportional to the square root of the diffusion coefficient, and the 
diffusion coefficient is proportional to E 2 when Larmor radii are larger than ~ L c [20]. We 
found no significant change to this conclusion when varying turbulent GMF parameters in 
the ranges tested in the previous section. However, the time spent in the source region for 
Bq = 8 /xG is ~ 30% smaller than for Bq = 4/xG at E/Z ~ 10 18 eV. For stronger turbulent 
fields, the larger turbulent field gradient towards z tends to make CRs leave the source region 
faster. 

We show in Fig. 7, which regions of the Galactic plane are passed most by CRs back- 
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traced from the Earth for E/Z = 10 18 eV/26 (upper left panel), 3 x 10 18 eV/26 (upper right), 
10 19 eV/26 (lower left) and 3 x 10 19 eV/26 (lower right). This equivalently shows which parts 
of the source region contribute most to the Galactic CR flux detected at Earth. The color 
code presents the fraction of all particles backtraced from the Earth which pass cubes of 
200 pc lateral size located in the Galactic disk. Here, CRs are not counted more than once. 
Location of Earth is marked by the bright spot at (x,y) = (0,8.5kpc). One can clearly 
see the shapes of the spiral arms present in the model of the regular GMF model. CRs in- 
deed diffuse or propagate faster along the regular field direction. The Galactic center region 
appears strongly demagnified, which makes it unlikely to significantly contribute to the ob- 
served fluxes. This is due to the stronger field in the disk and especially towards the Galactic 
bulge. For the forward tracking point of view, CRs potentially emitted by the Galactic center 
would escape from the Galactic thin disk which contains the Earth, and propagate towards 
larger z, before reaching Earth. At E/Z = 3 x 10 19 eV/26 (lower right panel), the CRs 
start to be in the ballistic regime and only the region within a few hundreds of parsecs from 
Earth could contain sources. At these rigidities, CRs cannot have Galactic origin any more 
because anisotropies would exceed the Pierre Auger upper limits on the dipole anisotropy, 
as discussed in the previous section. 

When the confinement time of CRs in the source region (~ 0.5 Myr for 1 EeV iron nuclei) 
starts to be comparable to the period between two potential Galactic UHECR sources, the 
continuous source distribution approximation used in this paper breaks down. For GRBs, 
this happens at E/Z > (0.1 — 1) EeV/26. In this case, the expected Galactic CR flux at Earth 
and its anisotropy would strongly vary on time scales of several Myr, see for instance Ref. [47]. 
The CR anisotropy at Earth may then be substantially smaller or larger than those computed 
for a continuous source distribution. Potential sources such as magnetars are expected to 
have a larger rate of ~ 1 per 1000 years. For such rates, ~ (0.5Myr)/(1000yr) ~ several 
hundreds of sources would contribute to the flux observed at E/Z = 10 18 eV/26. This larger 
number of sources substantially reduces the fluctuations in time of the Galactic CR flux and 
anisotropy at Earth. Therefore, in this case, the anisotropy does not significantly differ from 
the values presented in this work, and the continuous source distribution approximation is 
valid. 

Let us now estimate more quantitatively when the continuous source distribution ap- 
proximation breaks down. The average number of sources that would contribute to the 
Galactic CR flux at one given rigidity can be estimated as the average time spent by CRs in 
the source region multiplied by the source rate 1Z. We plot in Fig. 8 (left panel) this estimate 
of the number of contributing sources, for three different rates 1Z = 10, 1, 0.01 kyr -1 and 
for CR rigidities in the range E/Z = (1 — 32)EeV/26. We assume here for the turbulent 
GMF a Kolmogorov spectrum with L max = 200 pc, strength Bq = 4/iG and scale height 
zq = 2kpc. For sources with a rate 1Z = lkyr -1 (green line) similar to that expected for 
magnetars, the average number of contributing sources stays above > 100, up to the ankle 
for iron nuclei. It decreases from ~ 500 for 1 EeV iron to only ~ 10 at 32 EeV. As shown 
below, ~ 10 sources is too small for the continuous source distribution approximation to be 
valid. With a ten times larger source rate 7Z = 10 kyr -1 (red line), there would be > 100 
sources contributing for all the explored rigidity range, but 10kyr _1 is of the order of the 
Galactic supernovae rate which looks very unlikely for extreme CR accelerators. The blue 
line corresponds to rarer transients 1Z = 0.01 kyr" 1 . At the ankle, only one source would 
contribute in average. Sources such as Galactic GRBs with 1Z ~ 1 Myr -1 cannot be described 
by the continuous source distribution approximation. However, such sources are unlikely to 
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Figure 8. Left panel: Estimate for the number of sources that would contribute to the Galactic CR 
flux at Earth versus 26 x E/Z, for three different source rates 1Z = 1/(100 yr), l/(lkyr), 1/(100 kyr); 
Right panel: Estimate for the minimum rate TZ of Galactic CR sources that would be required to 
maintain relative fluctuations of the Galactic CR flux at Earth <j{F)/ (F) below re 5, 10, 25, 50, 100%, 
versus 26 x E/Z. For both panels, same Galactic magnetic field parameters as in Fig. 6. 



be responsible for the sub-ankle CR flux if it were to be of Galactic origin. Indeed, one can 
for example hardly match the bumpy CR spectrum resulting from rare Galactic transients 
to the observed smooth power law spectrum [47]. 

We provide in Fig. 8 (right panel) an estimate of the minimum rate 1Z of Galactic 
sources that would be required to maintain relative fluctuations of the flux cr(F)/ (F) below 
five given thresholds (~ 5, 10, 25, 50, 100%) on S> 1Z~ l time scales. Since we follow individual 
CR trajectories, we cannot directly provide ct(F)/ (F) for computing time reasons. We can 
however estimate its value: The panels of Fig. 7 also give an estimate of the total flux that 
would be received in any point of the Galactic disk from one source located at (x, y) = 
(0,8.5kpc). Looking for the total flux received at the Earth position from N sources -with 
same power- located in N random positions in the Galactic disk is then roughly equivalent 
to putting N observers in N random locations in the disk and summing up the total flux 
they receive from the single source located at (x,y) = (0, 8.5 kpc). We compute the flux 
F in the latter way for 10 4 different configurations, using the computations of Fig. 8 (left 
panel) for N. This yields the estimate for a(F)/ (F) that is used in Fig. 8 (right panel). 
To maintain <j(F)/ (F) below ~ 5% (red solid line) at the ankle for Galactic iron primaries, 
sources with rates comparable with that of Galactic supernovae 1Z ~ 30 kyr" 1 would be 
needed. For rates TZ ~ lkyr" 1 , a(F)/ (F) sa 10% for 1 EeV iron nuclei and remains < 25% 
below the ankle. In this case, the continuous source distribution approximation is still valid. 
For energies E > (10 — 20) EeV, it quickly starts to break down: For 20 EeV (resp. 30 EeV) 
iron nuclei, u(F)/ (F) 50% (resp. « 100%). For Galactic source rates 7Z < 0.01 kyr -1 , 
cr(F)/ (F) exceeds 100% above ~ 1 EeV and the anisotropy measured at Earth is expected 
to significantly differ from the averaged values presented in the previous section. 

6 Conclusions and Perspectives 

In this work we studied the consistency of a transition from Galactic to extragalactic CRs with 
existing anisotropy limits as a function of energy above E = 10 18 eV. The diffusion approx- 
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imation predicts a dipole anisotropy 5 = — 3Dj,-Vj ln(n) increasing with energy, since both 
the diffusion tensor Dij and the relative CR gradient Vj ln(ra) increase with energy. However, 
this approximation becomes unreliable at 0(E/Z) ~ 10 16 eV, and therefore we studied the 
propagation of CRs in the Galactic magnetic field directly by backtracking trajectories. We 
simulated the turbulent magnetic field on nested grids which allows one to include turbulent 
field modes B(k) with arbitrary small wave-lengths. For the regular Galactic magnetic field 
we used up-to-date models from Ref. [23]. Because the global structure of the GMF is still 
rather uncertain, we studied the dependence of the resulting anisotropy on the magnetic field 
parameters such as its strength Bq, scale height zq, correlation length L c and exponent a 
of its power-spectrum. We also examined the dependence of our results on the width and 
height of the disk in which sources are located. 

The main results of this study are presented in the Figs. 3-5. They show that the 
anisotropy mostly depends on the amplitude Bq of the magnetic field in the disk. As our main 
conclusion from this study, we found that existing anisotropy limits are not compatible with 
light (proton) and intermediate (CNO) nuclei of Galactic origin as dominant contribution to 
the CR flux above 1 EeV. By contrast, Galactic iron nuclei as CR primaries are consistent 
with the existing limits even up to 10-20 EeV, if the strength of the turbulent field is as large 
as B vms ~ 8 /xG. This finding implies that determining the chemical composition of the CR 
flux around 10 18 eV settles also the question of the transition energy between Galactic and 
extragalactic component: As light nuclei at this energy are not sufficiently isotropized, they 
have to be extragalactic. Therefore the fast increasing proton contribution indicated by the 
KASCADE-Grande collaboration between 10 eV and 10 18 eV suggests the beginning of an 
extragalactic component. 

We also studied qualitatively the dependence of the anisotropy on the effective density 
of sources, see Figs. 6-8. The average escape time of iron nuclei with 10 EeV energy from 
the Galaxy is ~ 10 5 yr. Assuming for magnetars a rate of 10 _3 /yr, the effective density 
of magnetars as sources of CR at 10 EeV is ~ 100/Galaxy. Thus magnetars satisfy the 
anisotropy constraint and can be natural candidates for the sources of the high-energy end 
of the Galactic CR flux in the scenario where the transition from Galactic to extragalactic 
cosmic rays occurs at the ankle, provided they are able to accelerate iron up to few x 10 18 eV. 

In summary, we conclude that models with a transition from Galactic to extragalactic 
cosmic rays around the ankle are consistent with the existing anisotropy limits if the compo- 
sition of Galactic cosmic rays at E > 10 18 eV is dominated by heavy nuclei. In contrast, if 
the chemical composition at these energies turns out to be light or intermediate, a transition 
at the ankle would be very strongly disfavoured. 
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A Validitation of the nested grid method 

In this appendix, we verify the validity of the new method we propose to generate turbulent 
magnetic fields on nested grids, see Section 2. We have reproduced the earlier results of 
Refs. [16, 20], and present below our computations for the diffusion coefficient versus the 
Casse et al. ones [20]. 

In Fig. 9, we present the numerical results for the parallel (right panel) and perpendic- 
ular (left panel) diffusion coefficients for 1 PeV to 600 PeV CR protons diffusing in a field 
containing both a regular and a turbulent component. Let us denote B reg the strength of the 
regular component, and B rms the root mean square strength of the turbulent one. Two levels 
of turbulence are used rj = 0.1 and rj = 0.46, with rj = B 2 ms / (B 2 cg + B 2 ms ). The magnetic 
field strength is set to 4 /j,G. For the turbulent component, we take a Kolmogorov spectrum 
(a = 5/3) with L max = 150 pc and L' min = 0.1 pc. Therefore, L' min is smaller than the Larmor 
radius tl for all energies. We average over a few turbulent field configurations and propagate 
2000 protons. 

Red symbols in Fig. 9 represent our results with the nested grid code introduced in 
Section 2, and green symbols represent the results of Figs. 4 and 5 of Ref. [20] adapted to 
the values we use here for L max and the magnetic field strength. The uncertainties of our 
values and those of Casse et al. [20] can be estimated from the fluctuations from one point 
to another compared to averaged behaviour of the diffusion coefficients. Both for rj = 0.1 
and for rj = 0.46, our results reproduce very well those of Casse et al.. We did not report in 
Fig. 9 the results of Casse et al. below E ~ 4 PeV because they correspond to computations 
for CRs with Larmor radius smaller than the minimum size of turbulent magnetic field 
fluctuations in their grid. This stresses one of the advantages of our new method: There 
is no lower limitation on L' min / X max because one can always add another grid with smaller 
spacing and smaller scales of the magnetic field fluctuations. Therefore one can safely explore 
low rigidities. 

In Fig. 10, we test our code for the case of 100 TeV to 1 EeV CR protons diffusing in a 
purely turbulent field, i.e. without any regular field, B Teg = 0. The parameters are the same 
as for Fig. 9, except for E = 100 — 300 TeV where we take L' min = 0.01 pc to ensure that L' min 
is smaller than 7*l- For computing time reasons, we keep L' min = 0.1 pc for E > 1 PeV. We 
use .B rms = 4/uG and a = 5/3. 

Red crosses in Fig. 10 correspond to our computations and green ones to those of Fig. 4 
of Casse et al. [20] adapted to our values for L max and B Tms . One can see that also in 
this case the results are in very good agreement. As predicted theoretically, the diffusion 
coefficient is proportional to E 1 ^ 3 (respectively to E 2 ) at low (respectively high) rigidities. 
The computations of Casse et al. for pure turbulence were done with a field generated 
superposing Fourier modes. This enabled them to check such a wide rigidity range. However, 
such computations are significantly slower than ours, as discussed in Section 2. 

For the computations at E > 1 PeV in Fig. 10, we take two grids and the intermediate 
scale L2 between the two grids is 5pc, which corresponds here to the Larmor radius of 
E ~ 10 PeV protons. We have computed twice more values in the range 3 — 30 PeV than 
at other energies so as to check that this scale does not induce any artificial imprint in the 
results. As can be seen in the figure, the diffusion coefficient behaviour in this energy range 
is smooth and not affected by that intermediate scale. 
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Figure 9. Computations of the perpendicular (left panel) and parallel (right panel) diffusion coeffi- 
cients for CR protons with energies E = f — 600 PeV, and for two different levels of turbulence 77 = 0.1 
and 77 = 0.46, a = 5/3 and L max = 150 pc for the turbulent magnetic field. 4^G for the magnetic 
field strength. Red symbols for our results and green ones for the Casse et al. results [20]. 
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Figure 10. Computations of the diffusion coefficient for CR protons with energies E = 100 TeV to 
1 EeV, in pure isotropic magnetic turbulence (no regular field); a = 5/3, B TIas = AfiG and L max = 
150 pc for the turbulent field parameters. Red crosses for our results and green ones for the Casse 
et al. results [20]. 
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